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DISCRETE-TIME MODELS FOR 
IMPLICIT PORT-HAMILTONIAN SYSTEMS 


FERNANDO CASTANOS*, HANNAH MICHALSKAt, DMITRY GROMOV*, AND VINCENT 

HAYWARD? 

Abstract. Implicit representations of finite-dimensional port-Hamiltonian systems are studied 
from the perspective of their use in numerical simulation and control design. Implicit representations 
arise when a system is modeled in Cartesian coordinates and when the system constraints are applied 
in the form of additional algebraic equations (the system model is in a differential-algebraic equation 
form). Such representations lend themselves better to sample-data approximations. Given an implicit 
representation of a port-Hamiltonian system it is shown how to construct a sampled-data model that 
preserves the port-Hamiltonian (PH) structure under sample and hold. 

Key words. port-Hamiltonian systems, nonlinear implicit systems, symplectic integration, 
sampled-data systems 


1. Introduction. The class of Hamiltonian systems has a prominent role in 
many disciplines. It was recently extended in [6] to include open systems, i.e. systems 
that interact with the environment via a set of inputs and outputs (called ports), 
giving rise to port-Hamiltonian (PH) systems. Such extended models immediately 
reveal the passive properties of the underlying systems, making them particularly 
well suited for designing passivity-based control (PBC) laws. Two types of model 
representations of Hamiltonian systems are in widespread use: the explicit repre¬ 
sentation stated in the form of an ordinary differential equation (ODE) on an ab¬ 
stract manifold [28{ l29l 125] 126] and the implicit representation stated in the form of a 
differential-algebraic equation (DAE) usually evolving in a Euclidean space [4]. While 
explicit representations are usually preferred in the context of analytical mechanics; 
see mmm, the implicit DAEs models lend themselves better for numerical compu¬ 
tations as they lead to simpler expressions for the Hamiltonian functions. The formal 
relations between the two representations and their equivalence can be established 
if the system’s configuration space is regarded as an embedded submanifold of the 
Euclidean space; see (4] for a full development. 

Of principal interest here will be the construction of sampled-data (discrete-time) 
models of PH systems that best approximate their continuous-time counterparts. 
Sampled-data models are important in digital control implementations and permit 
for simpler design of PBC laws directly in discrete time. In this context, the notion 
of “best approximation” deserves clarification. 

For linear systems, exact sampled-data models can be constructed by requiring 
that the solutions of the sampled-data and continuous-time systems coincide at the 
discretization points. Point-wise model matching is usually impossible in the case of 
nonlinear systems short of explicit derivation of analytical expressions for their solu¬ 
tions. For general dynamic systems, it is thus the choice of an integration method 
for generating an approximate solution (like Euler or Runge-Kutta) whose precision 
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relates to the compromise between complexity and order of approximation of the 
given continuous-time system by its sampled-data counterpart. In the special case of 
Hamiltonian systems, structure preservation is usually the main criterion for choosing 
a numerical method; see, e.g., m for more details on structure-preserving numerical 
schemes. Structure-preserving methods guarantee accuracy in long-horizon simula¬ 
tions. 

It is known that autonomous Hamiltonian systems conserve two quantities: the 
Hamiltonian function H (i.e., the energy or storage function) and a certain two- 
form ui, called the symplectic form. Numerical integration algorithms can either 
conserve H [Tans] or uj m, but not both. Conservation of u is usually preferred 
over conservation of H as the symplectic form unambiguously defines the class of 
Hamiltonian vector fields; see Theorem 12.31 and Remark |T| 

Contribution. Adopting a symplectic-form preserving approach, a sampled- 
data model for a given continuous-time PH system in developed here. The approach 
employs implicit representations of PH systems as the latter lead to simpler expres¬ 
sions for the Hamiltonian functions. Specifically, the Hamiltonian functions arising 
from implicit representations lend themselves well to the application of flow-splitting 
numerical integration methods; see |27] for a splitting method that applies to au¬ 
tonomous Hamiltonian systems (Hamiltonian systems without ports). Additionally, 
the discrete-time modeling approach presented here preserves passivity of PH systems 
lTheorem l4.2[) . The passive structure is preserved in the sense that the discrete model 
is the exact representation of another, possibly non affine, continuous-time PH sys¬ 
tem which, up to an approximation error of order two with respect to the sampling 
interval, has the same storage function H and the same output function y. 

Paper structure. Section [2] introduces the implicit representations of PH sys¬ 
tems. This section recalls the results from |4] and serves mainly to introduce the 
notation and to state the main assumptions ([T| and [2J. Also, the definition of PH sys¬ 
tems is extended to the case of non affine control systems. Splitting and symplectic 
integration methods for autonomous systems are recalled in § |3] Discrete-time mod¬ 
els for implicit PH systems based on vector splitting methods are developed in § [4] 
and their asymptotic properties are analyzed. An illustrative example is provided 
confirming the findings. 

2. Implicit representations of port-Hamiltonian systems. We begin this 
section by defining the configuration and phase spaces of Hamiltonian systems in 
implicit form. Conservative properties are also discussed. Similarly to the case of 
explicit representations of Hamiltonian systems, both the Hamiltonian function and 
the symplectic two-form of the implicit representation are invariant under the flow 
of the system. Following [6], input and output ports are added to the implicit rep¬ 
resentation, insuring passivity of the resulting PH system. Finally, we extend the 
definition of PH to the the case of non affine control systems. This will be needed 
when performing backward error analysis later in § [4] since time discretization entails 
the loss of affinity. 

2.1. The configuration space. The class of systems considered here is re¬ 
stricted to mechanical finite-dimensional systems with scleronomic constraints that 
evolve in continuous time. Systems of this type typically consist of M rigid bodies 
held together as one structure by the action of constraint forces. The position of each 
of the rigid bodies can be unambiguously described in terms of the position of its cen¬ 
ter of mass and its orientation in space. The configuration space of a spatial system as 
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a whole can thus be viewed as a subset, G, of a Cartesian space of dimension n = 6 M 
that is expressed in terms of smooth, independent constraint functions g : R™ —> R fc , 
with k < n, as the level set 

G := < 7 -1 ( 0 ) = {r£ R™ | g(r) = 0} (2.1) 

where g _1 (0) is the inverse image of 0 € R fc . Functional independence of the con¬ 
straints is expressed in terms of the following rank condition. 

Assumption 1. The rank of g is equal to k at all points of the set <? _1 (0). 

Recalling that the rank of a mapping is the rank of its tangent map (push- 
forward by g), then in any coordinates, the condition requires that the Jacobian 
G := {dg.jjdr 3 } !? has full rank, i.e., rank(G) = k at every point of <? _1 (0). The value 
0 G R fc is then called a regular value of g , the level set g _1 (0) is a regular level set of 
g , and g is said to be a defining map for g _1 (0); see [16], p. 113-114. It hence follows 
that G can be given the differentiable structure of a closed embedded submanifold of 
R”; see [16], Corollary 5.24, p. 114. The Constant-Rank Level Set Theorem (jl6j. 
Theorem 5.22, p. 113), further specifies its dimension as o = n — k. 

Because G is an embedded submanifold, it can also be regarded as an abstract 
o-dimensional manifold with local coordinates q (called generalized coordinates in 
mechanics). It is then also possible to exhibit an injective inclusion map: » : G 
R" (embedding), with rank(i) = dirnG, which is a homeomorphism onto the image 
*(G) C 1" and such that i(q) = r and g o % = 0. This inclusion serves as a map 
between local coordinates q and global Cartesian coordinates r. 

2.2. Hamiltonian equations. Let T* G be the cotangent bundle of G and let 
{llGPi} be local coordinates. A system is said to be Hamiltonian if its trajectories 
are integral curves of the Hamiltonian vector field D^ : T*G —r T(T* G), 

dH d dH d 

Fr dpi dq l dq' dpi 


where H : T *G —>• R is the sum of the system’s kinetic and potential energies K : 
T *G —> R and V : G —>• R, respectively; see [Tj for more details and a coordinate-free 
definition. 

The implicit model for a Hamiltonian system is defined as follows. Let T*R" 
be the cotangent bundle of R” and let { r l ,pi } be global coordinates. The implicit 
Hamiltonian vector field takes the form 


with 


A '.H.q = Dh — A 


dgi d 

dr 1 dpi 


9 = 0 , 


dH d dH d 
dpi dr F dr 1 dpi 


( 2 . 2 ) 


(2.3) 


the unconstrained part of the Hamiltonian vector field and H : T* R n —> R is again 
the sum of the kinetic and potential energies, but expressed in global coordinates. 
That is, K : T*R n -» R and V : 1" -> R; see pH 1211E] for details on the derivation 
of this equation. The vector field unravels as the DAE system 


r = V p H(r,p), p=-V r H(r,p) - G(r) r A 
0 = g(r) . 


(2.4a) 

(2.4b) 
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By applying Xfj g to both sides of the constraint equations g° = 0, one obtains 
the so-called hidden constraints 


f j ■■= X H , g (9 j ) 


dH dgi 
dpi dr i 


Thus, the system evolves on the closed submanifold 


(2.5) 


LG = {(r,p) € T*R n | g(r) = 0, f(r,p) = 0} (2.6) 

and we have Xfj g : LG —)• T(LG). A more rigorous construction of the phase space 
LG is given in [3], where it is regarded as an embedding of T*G in T*R". The map 
p i->- p as well as the formal relation between H and H can also be found in (4j . 

The Lagrange multipliers A j are defined implicitly by (12.211 . Precisely, application 
of Xn, g to the hidden constraints makes A j appear: 

= 0. (2.7) 

Thus, if the matrix 

f dg j df \ _ f dg j d 2 H dg l \ 

\ dr i dpi j jt { dr * dpidp m dr m j jt 

is non-singular on LG, then there are unique A j satisfying m and forcing the integral 
curves to stay on LG. 

Assumption 2. The Hessian matrix {d 2 H(r. p) / dpidpj } .. is positive definite for 
all ( r,p ) G T* R ra so H(r,p) is convex in p. 

This assumption, satisfied by most mechanical systems, together with Assump¬ 
tion [1] ensures that the matrix in (12.81) is invertible and the system is well defined. In 
mechanical systems, A is the covector of constraint forces that insure satisfaction of 
the constraints during the evolution of the system. 

2.3. Energy conservation and symplecticity. It is well known that the flows 
generated by Hamiltonian vector fields in explicit form preserve the Hamiltonian func¬ 
tion, i.e. the total energy of the system is conserved during the evolution of the sys¬ 
tem [la m Ei]. As the explicit and implicit system representations are equivalent, 
the conservation also holds for the implicit Hamiltonian vector field (12.21) , as is easily 
confirmed by computing the Lie derivative of H along the flow generated by Xn, g HD- 
Indeed, 

Cx a ,H = X H , g (H) = D h (H) A^g = A j f j , (2.9) 

where the third equality holds because Dh{H) = 0 and because of the definition of 
/ given in (12.51) . Eq. (12.91) shows that Cx H<g H | iG = 0 (cf. Eq. (12.61) 1. so H remains 
constant along the system trajectories. 

It is also well known that, besides the 0-form H 1 Hamiltonian flows preserve a 
certain 2-form called the symplectic form. In the case of implicit representations the 
symplectic form is given by the formuleS 

w := dr* A d pi , (2-10) 


1 See [T] for a coordinate-free definition. 
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which acts on vectors of T(T*R”), with Einstein’s summation convention implied. 
Definition 2.1. A differentiable mapping f : T*R" —» T*R" is called symplectic 

if 

f*U}=U. ( 2 . 11 ) 

Here, </>* denotes the pull-back map associated with f, defined by 

Definition 2.2. Let <f> : Mi —> M .2 be a smooth map of manifolds and let 
p G Mi. The pull-back map f* : TL^AA 2 —> T* M 2 associated with <f> is a dual map 
to the push-forward map <j>* and is characterized by 

(<!>*€,x) = {fff*x) for f G x g t p m 1 . 

The application of this definition permits to re-write (12.111) in the equivalent form 

w(</>*£, 0 * 77 ) = w(£, 77 ) for all £,77 G T(T*R") . (2-12) 

The conservation of w along Xn, g can be established by showing that the Lie 
derivative Cx H 3 w is equal to zero, the demonstration bearing similarity to that of the 
conservation of H; see [5] for the explicit derivation. Here, we cite the result 

Theorem 2.3. mm Let H be twice continuously differentiable. The flow 
(ft ■ LG —> LG of Xn, g governed by m is a symplectic transformation on LG, i.e., 

(fffuJ = LO 


for every t for which ft is defined. 

Remark 1 . The converse statement, that every symplectic flow <j> t solves Hamil¬ 
ton’s equations for some H, is also true, so symplecticity is a characteristic property 
of Hamiltonian systems W- This does not translate to the case of energy conservation, 
i.e., while every Hamiltonian system conserves energy, not every energy-conserving 
system is Hamiltonian. 

2.4. port-Hamiltonian systems. In the presence of external forces and dissi¬ 
pation it is convenient to represent ( 1 ^ 1 ) as an input-output system equipped with a 
pair of port variables (u, y), giving rise to a PH system ; see [2D11521125] for the original 
definition as stated with respect to Hamiltonian systems in explicit form. Extend¬ 
ing on this definition the Hamiltonian systems in implicit form, a port-Hamiltonian 
system is described in terms of the vector field XH,u,g '■ LG x (R m )* —>• T(LG): 


X}J,u,g D 


H 


(uiUf - Xj 


d_l_ 

dr 1 


_d_ 

dpi 


9 = 0 


(2.13) 


with u G (R m )* defined as the controlled or input variable, y G R m is defined as the 
output variable that satisfies 


y‘ = u‘ a X 

dpi 


and where Uf are maps from R” to R. 

The vector field (12.131) and the output unravel as the DAE system 


f = V p H(r,p) , p = —X r H(r,p) — G(r) T A + U{r)u 
y = U(r) T X7 p H(r,p) 

0 = g(r) . 


(2.14) 
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By analogy with the results described in Sec. 12.21 one can determine the Lagrange 
multipliers A explicitly. The constraints f a = 0 imply that 

x H,u, g (f a ) = D H (f a ) + = (2.15) 

from where it follows that, as long as (12.81) is non-singular, there are unique A j (in 
general dependent on u as well as on r and p ) such that XH tU ,g(f a ) = 0 and such that 
the integral curve stays on LG. 

It can be readily seen that an implicit PH system described by (12.131) no longer 
preserves H. The Lie derivative of H is now 

£x h , u JH) = X H , u , g (H) = DnW+mU, 1 — - A^— = w - Xj f j , 

which establishes the power balance 

Cx H , u , g H\ LG = u iy l . (2.16) 

Since the product uiy 1 is equal to the rate of change in energy, we say that (it, y ) is a 
power-conjugated pair of port variables. If, in addition, the restriction of H to LG is 
bounded from below, i.e., if the image of LG under H is bounded from below, then 
(12.131) is called passive, or more precisely, lossless. Boundedness of H can be easily 
assessed using the following proposition. 

Proposition 2.4. $ If the potential energy V is lower semi-continuous and G is 
compact, then H restricted to LG is bounded from below (hence, the vector field (12.131) 
describes a lossless system). 

With the inclusion of the control variable u, it can no longer be expected that 
the flow of (12.131) be symplectic. Indeed, it is not hard to see that the Lie derivative 
of u> along X}j,u,g is in general different from zero. 

Proposition 2.5. m The Lie derivative of u> (12.101) along (12.131) and restricted 
to LG satisfies 

£x H , u ,M LG =dr i /\d(u l U i l ). (2.17) 


Example: A double planar pendulum. Let us recall an example from [4], on 
which we will elaborate when discussing sampled-data models. 

Consider the model of a double planar pendulum shown in Fig. l2.1l that comprises 
a pair of point masses m a and mb whose coordinate positions are r a = ( r ax , r ay ) and 
r b = (r bx ,r by ), respectively. The massless bars are of fixed lengths l a and lb which 
gives rise to the two holonomic constraints: 

<?>) = ||r°|| 2 — l 2 a = 0, g\r) = \\r s f-ll = d (2.18) 

with r := ( r a ,r b ) G R n , n = 4, k = 2, r s := r b — r a . The rank of the constraint 
Jacobian is full as 

/ JIG X rpCLy Q 0 \ 

rank G(r) = rank _ f Sy ^ r s y J=k (2.19) 

for all r G G. Therefore, 0 is a regular value of g and G is an embedded submanifold 
of R 4 . 
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Fig. 2.1. A double planar pendulum, a simple Hamiltonian system. 


The total energy is 


H(r,p) = -p T M 1 p + g(m a r ay + m a r ay ) , 


where M := 


'rn a I n On 

On nTfolyj 

Substituting (12.181) and (12.201) in (12.41) gives 


and g is the standard gravity. 


r a = m a 1 p a , r b = m b 1 p b 


(Pa x \ 


0 j 


/ r ax 

- rpbx \ 

Pay 

Pb x 

= 

gm a 

0 

-2 

y*a y 

0 

— jAy 

\Pby) 


\gm b ) 


\o 

r Sy J 


( 2 . 20 ) 


(2.21a) 

(2.21b) 


which, together with (12.181) . constitutes a set of DAEs describing the motion of the 
double pendulum in implicit form. The multipliers Ai and A 2 are the magnitudes of 
the tension along the two bars. 

Now, assume that the double pendulum is actuated by application of torques 
rti and U 2 to the joints that correspond to the angles q 1 and q 2 , respectively. The 
resulting linear forces are then U 1 iti and U 2 U2 with 


U 1 :={U i 1 } i = 


(—r‘ 

0 
0 


1 


/ r 5y \ 


and U 2 := {U 2 }, = 


V 0 / 


—r 
A 


\ r° x ) 


— U l . 


( 2 . 22 ) 


The manifold defined by (12.181) is compact and the potential energy is continu¬ 
ous, which confirms that the double pendulum is passive with passive outputs y l = 

77 l _ tt l-i 

U i dpi — u i 1 ■ 

An explicit model for the double pendulum as an ODE can be derived by choosing 
the generalized coordinates on G as q 1 € (—7r,7r) and q 2 G (—7r, 7 t), as motivated by 
Fig. 12.11 The associated embedding r = i(q) satisfying g o 1 = 0 is readily exhibited 
as 


/ y*a x ^ 


/ l a cos q 1 \ 

f a y 


la sing 1 

y*bx 


l a cos q 1 + Zbcosg* 

\r by ) 


\l a sing 1 + l b sing 1 J 


(2.23) 
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In local coordinates, the total energy is 


H{q, p) = -p T M(q) 1 p + g ( m t l a sin q 1 + m b l b sin g‘) . 


with 

mi,) = ( m *'° + ,, 2 

v m b l(; + m b l a l b cos q z 
The motion of the system is described by 


nrr ^ + 2 mblJb cos q 2 m b ll + m b l a l b cos q‘ 

M \q) - 1 2 7 , m& ;2 


q = M(q) 1 p , p = -W q V{q) - V g ( )-p T M(q) 1 p)+u 


(2.24) 


(2.25) 


(see [I] for more details). 

Two representations for the same system were derived, one in the form of an 
ODE (12.2511 and the other as a DAE (12.2111 . The main point of this example can be 
summarized in the following remark. 

Remark 2. The implicit DAE representation of the Hamiltonian systems renders 
a separable Hamiltonian function (12.201) . i.e. the kinetic energy does not depend on 
r and the inertia matrix M is a constant diagonal matrix. Moreover, the potential 
energy is linear, which results in a constant potential energy gradient. On the other 
hand it is easily verified that the explicit Hamiltonian representation does not share 
such fortunate properties. The explicit Hamiltonian is no longer separable and the 
inertia matrix depends on the generalized coordinates q. 

The a cost of a simpler Hamiltonian function is the higher dimensional implicit 
model representation and the appearance of the Lagrange multipliers. As will soon be 
seen, however, implicit representations are particularly advantageous for the purposes 
of system discretization. 

2.5. Non affine port-Hamiltonian systems, ft was assumed in (12.1011 that 
the control variables enter the vector field affinely. Many physical systems exhibit 
this property so, from a modeling point of view, this assumption is not too restrictive. 
However, for the purpose of performing backward error analysis in §01 we will need to 
consider PH systems for which the control might enter in a non affine way. Motivated 
by the fact that £x HtU:g w\ LG = 0 when u = 0, we propose the following extended 
definition of a PH system. 

Definition 2.6. A controlled vector field X : LG x (R m )* —> T(LG) is said to 
be port-Hamiltonian if u = 0 implies that the generated flow is symplectic. 

Passivity of a non affine port-Hamiltonian system can be established by redefining 
the passive output. 

Lemma 2.7. A (not necessarily affine) smooth PH system described by a vector 
field X can always be decomposed as X = Xo + uiZ 1 , where Xq : LG —> T(LG) is a 
Hamiltonian vector field and Z l : LG x (R m )* —>- T(LG) are the input vector fields. 
Hence, X satisfies the power balance X(H) = uiy 1 for some real-valued function H 
and real-valued output functions y l = Z l (H). (The output functions may now depend 
directly on u as well as on r and p.) 

Proof. Following [18], we first show that a smooth control vector field can be split 
into a drift and a set of vector fields having u factored out. Let us define the drift Xq 
as X 0 (x) = X(x,0) and let us define the vector fieids W l by the equations 

W\a) = -^-(X(a)) for all a e C°°(LG,K) . 
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It follows from the chain rule that 



Upon integration on both sides of the equation we arrive at 


L 



0 



(2.26) 



W l (x, 6u)(a)d6 for all a G C°°(LG,]R) . 


It follows from the hypothesis that Xo generates a symplectic flow, so it is a 
Hamiltonian vector field and satisfies Xq(H) = 0 for some real-valued function H. 
Applying X to H shows that X(H) = Xq(H) + uiZ l (H) = uiy 1 with y l = Z l (H). □ 
Remark 3. For an affine PH system, the formulae of this lemma recover the 
output functions (12.1411 with Z l = Ujthat is, y l = Z l {H) = 

3. Sampled-data models for autonomous Hamiltonian systems. Com¬ 
puting a sampled-data model of a dynamical system basically amounts to computing 
an approximate solution of the differential equations during a small interval of time. 
This problem has been studied extensively in the literature of numerical analysis, from 
which we borrow some results and terminology. In numerical analysis, a sampled-data 
model is the central component of an integration method or a numerical integrator , 
up to the point that these terms are used interchangeably. 

Mathematical models for sampled-data systems arise in diverse circumstances. In 
the direct approach to digital control, i.e., as opposed to the emulation of continuous 
control laws, the design of the controller is performed in discrete time, the designer 
working directly over a sampled-data model. When designing directly in discrete time, 
the controller can be directly implemented on a digital device. Also, it is possible to 
exploit the advantages of switched controls or, e.g., multirate control techniques m- 
Computing the sampled-data model for a given nonlinear system relies on the 
computation of a solution <\> t of the corresponding ODE or DAE, which is in general 
impossible to do analytically, so one has to settle for an approximate solution. 

When simulating the behavior of dynamic systems, a discrete-time model of the 
continuous system is also used for computing a numerical solution to the initial-value 
problem. Many different integration methods (or methods for short) can be found in 
the literature. Let us first focus on integration methods for autonomous systems and 
further restrict our attention to one-step methods defined by a transformation 


i/j h :x a ^- x a+ i , 

where the constant step-size h is regarded as a parameter of the methoc0. For a given 


initial condition xq in the phase space, ^h is applied recursively to generate a discrete 


2 For a more general method, the value of the x^+x need not depend only on x a , but may also 
depend on the previous values x a — 1 , x a — 2 ,- • • (a multistep method). Also, the value of h need not 
be constant in general. 
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flow xi,X 2 ,X 3 ,... that approximates the true flow <ph{xo), <p 2 h( x o), cp 3 h( x o), ■ ■ ■ of 
a given vector field X at time instants h, 2 h, 3 h, ... In this sense, the map tp^ is a 
discrete-time approximation of <ph (or a sampled-data model of (ph). 

Definition 3.1. A one-step method has order s if the local error satisfied 

tph( x o) - <ph{xo) = 0{h s+1 ) as h —)■ 0 (3.1) 

uniformly in x o. A one-step method is said to be consistent if s > 1. 

Let us now discuss some important properties of numerical integrators, like order 
and symmetry. 

3.1. Symplectic methods. If a sampled-data model approximates the discrete 
time behavior of a Hamiltonian system, one could hope for tph to inherit its funda¬ 
mental qualitative properties: energy conservation and symplecticity. Unfortunately, 
it is not possible to preserve H and to simultaneously, unless iph agrees with the exact 
flow cf>h up to a reparametrization of time [7). For this reason, one has to choose ei¬ 
ther in favor of one or the other invariant^]. Energy conserving methods have received 
some attention [9] |3T] [TO] [12] 13 , but in light of Remark [Q most of the literature 
focuses on symplectic integration algorithms (see pTJ [IT} [5] and references therein). 
A comparison between both approaches is carried out in m for the rigid body. 

A theoretical advantage of constructing a symplectic one-step method is that, 
even though iph only approximates (ph up to the s’th order, it coincides exactly (if 
one disregards convergence issues) with the flow of another Hamiltonian system, a 
modified Hamiltonian system described by a modified differential equation. 

Theorem 3.2. fill p. 352] A symplectic method iph : LG —>• LG for the con¬ 
strained Hamiltonian system (El) has a modified equation that is locally of the form 


r =+'V p H(r,p) , p = -V r H(r,p) - G(r) T A (3.2a) 

0 = g(r) (3.2b) 


with H = H + hH 2 + h 2 Hz + ... Furthermore, 


dHj (r , p) dg l (r) 
dpi dr i 


for all ( r,p) £ LG , 


all l = 1 and all j. Note that the actual value of the Legendre multipliers A 

differ in general from, those obtained for the original system 1 2-4[ ). 

In other words, for an initial condition xo, iph(xo) is equal to the solution of (13.21) 
at time t = h. Note that m provides information about the difference between 
the actual flow (ph of X and the approximate discrete flow iph. This is the kind 
of information that forward error analysis aims at. While certainly useful as an 
indicator of the quality of the approximation, Eq. m only evaluates the behavior 
of the approximate flow on the first iteration, but says nothing about its long time 
behavior. From (ED alone we cannot infer anything about the error x a — (j) a h( x o) 
when a is large, so we do not know if errors accumulate or if they average out to zero. 

3 We use big-O notation when quantifying approximation errors, i.e., for a given pair of functions 
e 2 (h), we write e\(h) = 0 (e 2 (h)) as h —>■ a as shorthand for limsup^v a || ei ^j|| < oo. 

ll e 2v.'OII 

4 For particular Hamiltonians there might be other invariants, such as momentum or angular 
momentum, but in general there need not be. 
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On the other hand, Theorem 1531 tells us that if 0 ^ is symplectic, then there 
exists a modified continuous system whose flow coincides exactly with the discrete 
flow generated by 0ft. The modified system ( 13 . 21 ) preserves the Hamiltonian structure 
of the original system ( 12 . 41 ) and it is ‘close’ to it in the sense that H = H + <D(h s ) 
for a method of order s. In other words, a symplectic integration method preserves 
the original 2-form oj and a different (but close) Hamiltonian function. This property 
guaranties that the good behavior of the integration scheme is maintained during 
many iterations, giving a global nature to the local property m■ This observation 
is at the center of backward error analysis m- 

3.2. Splitting methods. A practical advantage of symplectic schemes is that 
they lend themselves well to the application of splitting methods. To illustrate the 
idea, consider again the unconstrained or, otherwise, explicit Hamiltonian vector field 
D^ on an abstract manifold T*G. If the Hamiltonian function is separable, i.e., if it 
can be written as H(q,p) = H a (q) + Hb{p), then the vector field can be spilt into two 
Hamiltonian vector fields 

__dthd_ _dH h d 

dq l dpi a dpi dq l 

with 0^ = 0^ +Dfi b - Notice that, taken separately, each vector field can be trivially 
integrated exactly. For ( q a ,Pa ) € T*G (we use the subindex a to refer to an element 
in a sequence, not a particular coordinate) we have 

Q = (~V q H a (q)) Ow) = (pa - h • y q H a (q a )) = ^ h (X) 

and 

q a T h • ^7pHb(p a )\ 0 ^ j 

Pa ) ~ h ' h 

A first-order symplectic method can be easily constructed by performing the compo¬ 
sition 


(0\ _ f +Vp.Hb(p)'\ fda+l \ _ 

\PJ V o ) \p a+1 ) - 



— 0b ,h ° 0a ,h • (^-^) 

Indeed, the maps 0b,ft, and 0 a> ft are symplectic because they are the exact flows of 
Hamiltonian vector fields. Since the composition of two symplectic maps is again 
symplectic, 0ft is symplectic. 

Many simple Hamiltonian systems with phase space T*G are not governed by 
separable Hamiltonians, so splitting methods cannot be applied directly. However, 
the Hamiltonian function of many mechanical systems becomes separable if the phase 
space is embedded in T*R” (see, e.g., Remark [2]). An interesting symplectic method 
that is particularly well suited for this class of systems was proposed in Ezj. Roughly 
speaking, the idea is to compute a symplectic method for the unconstrained Hamilto¬ 
nian vector field Dh '■ T*M n —>• T(T*R ra ), i.e., a symplectic map ipH,h approximating 
the solution of the ODE (notice the absence of the constraint equations) 


at time t = h. 


r = +V p if(r,p) 
V = - V r H{r,p) 
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If H is separable, ipu,h can be readily found. The method T H,g,h for the original 
constrained Hamiltonian vector field Xjj g is then constructed by taking the image of 
ipH,h and applying a correction term that ensures that the value of tk H,g,h belongs to 
LG, so that the constraints are satisfied. The correction is done in a careful way so 
that the resulting map is still symplectic (see also HZ!)- Depending on the accuracy 
of tpH,h, the resulting i &H,g,h can be of first or second order (see § [4] for details). 

3.3. Symmetric methods. For each t for which the solution is defined, the 
flow <f>t(x o) of an autonomous differential equation defines a transformation on the 
phase space. It follows from the group property of the flow [2] that the inverse of 
the transformation can be obtained simply by reversing time, that is, (pf 1 {x\ ) = 
4>-t(xi) = xo■ Needless to say, this property does not hold in general for a discrete 
approximation iph, which motivates the following definition. 

Definition 3.3. The adjoint method ip j* of a method iph is the inverse map of 
the original method with reversed time step —h, i.e., 

r h := {ip-h)- 1 . 

In other words, x\ = ip^{x o) is implicitly defined by ip-h(x\) = xq■ A method for 
which ipfo = iph is called symmetric. 

From a theoretical point of view, an approximate discrete-time flow should be 
symmetric because actual continuous flows are. But symmetry is important from a 
practical point of view too. It has been proved in (54] that all symmetric methods 
are of even order, a fact that can be exploited to construct high-order methods from 
simple lower-order methods. For example, one can take a first-order non-symmetric 
method, compute its adjoint and construct a symmetric method 

’F/i = iph. o ip* h . (3.4) 

2 2 

We know that H/h is at least first order, but since 'i’h is symmetric, we also know that 
the order has to be even, so we conclude that the method is actually of second order. 

The scheme (13.411 works particularly well for splitting methods. Take, e.g., the 
integration scheme m - The maps (pb,h and (p & ,h are symmetric (because they are 
exact solutions of a differential equation), but their composition is not symmetric in 
general. To remedy this, one can compute the adjoint method 

'0/i = (0b,— h 0 <Pa., — h) = ( P a ^ — h ° < / > b , — h = ° 4>h,h (3-5) 

and, using (13731) . (1331) and m , construct 


^ h — (pb, | 0 ,! 0 ° — 0b,| 0 0a ,h 0 0b,| > 


which is a second-order symmetric method. 

3.4. Modified vector fields and exponential representations. Consider a 
vector-valued function F and a vector field X, both defined on LG. If F and A' are 
analytic, then the composition of F and the generated flow (pt(x o) can be expanded 
in a Taylor series around t = 0, 


00 ** 

F o (p t {x 0 ) = exp(tX)F{x 0 ) ■■= -j Xl ( F ){ x o) , 

i=0 *' 
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Fig. 4.1. Sampled-data PH system with sampling period h. The zero-order hold produces piece- 
wise constant inputs ui(t) = u la for t E [ah, ah-\-h), which is fed to the continuous-time PH system. 
The output is then sampled to generate the discrete-time output sequences {y l a }- 


where X°(F ) = F, X 2 (F) = X{X{F)), X 3 (F) = X(X 2 (F)) etc (see (22] |33] for 
details). In particular, if F is taken as the identity function Id, one obtains the flow 
(j)t{x o) = exp(fX) Id(xo). Since an s-order method i/'ft for X coincides with the flow of 
a modified vector field X = X + 0(h s ) JTTj, p. 340], it is also possible to expand il>h in 
a Taylor series, tph{x o) = exp(hX ) Id(xo). This exponential notation is a convenient 
way to express the relationship between a vector field and the flow generated by it, 
as well as to analyze the composition of flows. 

4. A splitting method for implicit port-Hamiltonian systems. Suppose 
that there is a sequence of commands {iq Q } Q gN. Each command in the sequence ar¬ 
rives at the discrete instants of time a = 0, h,2h,... where h is a positive real number 
—such commands could be generated, e.g., by a computer program. Suppose further 
that a zero-order hold transforms this sequence into piece-wise constant controls 

ui{t) = u la for t £ [ah, ah + h) , (4-1) 

which are fed into a PH system. Let cj)t(xo,u( ■)) be the integral curve of the non- 
autonomous vector field (12.1311 with control (14.11) and passing through xo £ LG at 
t = 0. Let 

y l {t) = °M x o,u(-)) 

be the corresponding outputs and let {y l a } with y l a := y l (ah) be the sequence ob¬ 
tained by sampling them at discrete instants of time ah (see Fig. 14.11) . We call the 
resulting system a sampled-data port-Hamiltonian system. 

The goal here is to develop a method for derivation of discrete-time (or sampled- 
data) models for PH systems given by implicit vector fields. The underlying idea 
is to split the PH vector field into two components: the vector field describing an 
unconstrained system with state-space equal to the whole T* R", and a vector field 
containing the Lagrange multipliers, the one that maintains the trajectories on the 
submanifold LG. Splitting the vector field simplifies the computation of the sampled- 
data model by decomposing the problem into two simpler subproblems. 

Now we extend the results of m to the PH case. We show that, with a straightfor¬ 
ward modification, the method presented in originally intended as an integration 
scheme for autonomous Hamiltonian systems, can be used to compute sampled-data 
models that preserve the main properties of a PH system. 

Consider again the implicit vector field 
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defined on LG, with Dh as in (12.31) and with piece-wise constant controls (14.11) . Sup¬ 
pose that a method ipH,u,h '■ T* R" —>• T* R ra of order s > 1 for the unconstrained 
vector field Xh, u = Dh + u iUi l -^p has been computed. Again, in many cases H is 
separable so a high-order and symmetric method with a PH modified vector held can 
be easily found. The controls ui are constant during each sampling interval, which 
further simplifies the task of Ending ipH,u,h- 


Let us dehne the map 


(Z) : = L - ^) t a) - (4 - 3) 

Loosely, this is an approximation of the integral curve of the remnant vector held 
—Aevaluated at t = h and subject to g = 0. More precisely, for arbitrary 
functions A ? of r and p. we have that 

D djXjgi) d d(Xjg J ) d -dA,- d f dgp i dx A d 

Xi93 dpt dr 1 dri dp 1 dpi dr 1 \ 3 dn dn) dp 1 

and, for r £ G, the vector held reduces to 


D x 


•i 9- 


= XD ■ = _ A .^1 
3 93 dr 1 dpi 


(4.4) 


In other words, when restricted to G, the vector held — A j OjL is Hamiltonian (hence 

it generates a symplectic how). 

Lemma 4.1. \27 | / Let g(r a ) = 0. Then, the map (14.31) is a first-order symplectic 
method for D^. g j. That is, for r a € G, 


$A ,h 


= exp(hD\ g o) Id 


D\ jg i — D 


A i9 i 


with A a modified or perturbed version of A. 

A method for (14.21) can be obtained from the symmetric composition 


^H,u,g,h = $^1 O l/) H ,u,h 0 \ ■ (4.5) 

For each ( r a ,p a ) € LG, the values of p and v are determined implicitly by the 
constraints g(r a +±) = 0 and f(r a +i,Pa+i) = 0 (i.e., by (r a +i,p a +i) £ LG). In this 
way, T H,u.g,h dehnes a transformation on LG. 

The transformation '&H,u,g,h produces an approximate discrete how for a given 
command sequence {u la }. From this how, an approximate output sequence can be 
obtained by evaluating the output function y l = at each discrete time ah. 

Theorem 4.2. Consider the implicit method '&H,u,g,h 63 and let Xh, u be the 
modified vector field ofipH,u,h- 

(i) The method preserves the constraints g J = 0, f J = 0 and is of order s = 
min(s, 2), where s is the order of 4>H,u,h- 

(ii) The method is symmetric if ijjH,u,h is symmetric. 

(iii) IfipH,u,h Is symplectic for u la = 0 (i.e., if Xh, u is port-Hamiltonian), then 
the modified vector field Xn,u, g ■ LG —> T(LG) is also port-Hamiltonian with Hamil¬ 
tonian and output functions 

H = H + 0{h s ) and y l = y l + Q(h B ) . 


(4.6) 
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Proof. The method preserves the constraints by construction. The proof about 
the order of the method follows the same lines as the one given in E3 except that, 
since we are dealing with port-Hamiltonian vector fields, Lie brackets have to be 
used instead of Poisson brackets. We will compute Xn, u ,g, the modified vector field 
generating 'F H,u,g,h , and show that it agrees with Xn,u,g up to the first or second 
order, depending on whether 'F H,u,g,h is, respectively, first or second order. 

Let us consider the case s = 1. Using the exponential notation and Lemma 14.11 
the composition (14.51) takes the form 


’F H,u, g ,h = exp 



Id, 


where Xh,u = Xh, u + 0(h). Applying the Baker-Campbell-Hausdorff (BCH) for¬ 
mula to the product of the first two factors and truncating after the first term 
gives 

exp (^D 0 . g ^j exp (hX H ,u^j = exp (hX'^j (4.7) 

with X' = Xh,u + + 0(h). Applying BCH again to include the third factor 

gives 


exp (hX'^j exp = exp ( 'hX H ,u , g ) (4.8) 

with the modified vector field Xn, u ,g = Xh, u + \ (D 0 . g j + D^. g j} +0(h). Using (14.41) 
and Xh,u = Xjj, u + O(h), we can write the modified vector field as 

Xu,u,g = X H , U + D g j + 0(h) . (4-9) 

The hidden constraints f l = 0 imply that 

X H ,uM l ) = XnAf) + if 1 ) + 0(h) = 0 . (4.10) 

It follows from (14.101) . (14.41) and (12.151) . that the Lagrange multipliers A, and the 
‘modified Lagrange multipliers’ i/j and fij are related by the equation (vj + fij )/2 = 
A j + 0{h), which when substituted back in (14.91) gives the desired result: 

X-H,u,g = Xh,u + ^jDgj + 0(h) = Xjj,u,g + 0{h) . 


For s = 2 we follow the same procedure, but we truncate the BCH formula after 
the second term. For the expression gZZD , the intermediate vector field is 


X’ = X 


H,u 


-D 


Vj9 j 


h - 
4 . 


D 


C'j g j > 


Xh,v 


+ (D(h 2 ) 


where [ •, • ] is the standard Lie bracket. Using the initial assumption Xh, u = Xu u + 
0(h 2 ), we can write X' as 

X' = X H:U + — Dp jg j + — \D 0 . g j,X H ^y\ + 0(h 2 ) . 
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Regarding (|4.8|) , the modified vector field for the complete scheme is 


Xh.u.o — Xh.i 


( D 


Vj9 


D 


H 9 J 


) 

h r -| h 

+ 4 [ D P jg J’ X H,u] + J 


Xh.i 


2^ C 'j9 j ’ 


0(h 2 ) . 


Using (14.41) and the skew symmetry and bilinearity of the Lie bracket, the vector field 
can be equivalently written as 


X-H,u,g = X H ,u + ~ 3 + ^ J Dgj + — [( Vj - fij)Dgj , Xh,u\ 

+ g \vjDgJ ,p>jD g j~\ + 0(h 2 ) . (4.11) 

In order to extract information from the equation Xjj iUt g(gj = 0, let us first open 
the brackets in (14.111) and write 


~ zb-+/L- h 

Xjl,U,g XfJ^ U ~f" “ Dgj — { Vj [ij)Dgj Xf{ u 

— ~Xh,U {{pj — L l j)Dgj) + — VjDgj (fijDgJ ) — (VjDgj^ O ( k ) . 

Taking into account that f l = XH,u(g l ) = 0 and D g j(g l ) = 0, we have that 

X H ,u,g(9 l ) = ~ N) D gif l + 0(h 2 ) = 0 , 

from which we can see that modified Lagrange multipliers satisfy the order relation 

D j -fi j = 0(h). (4.12) 

By substituting (14.121) back in (14.111) we can verify that the commutators are actually 
second order, that is, 

XH,u,g = Xh, u + 3 Dgj + — \pjDgj , (Vj + 0(h )) D g j^ + 0(h 2 ) 

Z o 

= X HtU + + op 2 ) . (4.13) 

From XH,u, g (f l ) = 0 and (12.151) we conclude that, when s = 2, (i>j + jij)/ 2 = A j + 
0(h 2 ), so the desired result follows: X HtU>g = X Hu +\jD g j +0(h 2 ) = X H ^ g +0(h 2 ). 

For statement (ii), notice that, when restricted to LG, the method (1X51) can be 
described by the implicit equations 

r a +l \ , ( r a 

p<x +1 + hG(r a+1 ) T p) H,u,h \p a - hG(r a ) T v 

g{r a +i) = g(r a ) 

f{r a +i,p a +i) = f(r a ,p a ) , 

where r a ,p a are the independent variables and r a +nPa+i are the dependent variables. 
The vectors v, p are (also dependent) dummy variables that can be discarded after 
r a +i,p a +i have been found. 


(4.14a) 

(4.14b) 

(4.14c) 
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After reversing time (that is, after substituting h by —h), equation (14.14al) be¬ 
comes 


r a +i 

P a +1 - hG(r a+1 ) T p 


= H,u,-h 


p a + hG(r a ) T v 


Recall that ipH,u,-h = u, h ^ ’’Ph.u.ii is symmetric. Therefore, when restricted to 
LG. the reverse-time method is 


Pa + hG(r a ) T v) 

g(r a +i) = g{r a ) 
f{r a+1 ,p a +i) = f(r a ,p a ) 


< a+1 


p a +1 - hG[r a + 1) /i 


(4.15a) 

(4.15b) 

(4.15c) 


which is the same as (14.141) . but with r a ,p a and v interchanged with r a+ i,p a+ i and 
p, respectively. This implies that, if we input r a +i,p a +i as independent variables, we 
recover r a ,p a as the dependent variables, that is: 4' H,u,g,-h is the inverse mapping of 
^H.u.g.h- (In general, the vectors v , p obtained using (14.141) will be different from those 
obtained using (14.151) . but this is inconsequential since these are dummy variables.) 

In statement (iii), the fact that XH.u.g is PH follows directly from Lemma 14.11 
Definition 12.61 and the fact that the composition of symplectic maps is again symplec- 
tic. In other words, '&H,u,g,h is symplectic when u la = 0, so 


X H.u.g — Xg + UlZ 1 — Dfj + UlZ 1 + y Dgj 


(4.16) 


for some Hamiltonian function H and some input vector fields Z l . Since the method 
is of order s, we have 


Xh.u.q = Dh 


u l U l l — + \W gj +0(h s ) . 
dpi 


(4.17) 


By setting ui = 0 and recalling that A j = A j+0(h s ). it follows from (14.161) and (14.171) 
that Dfi = Dh + 0(h s ). which implies that H = H + (D(h s ). and, in turn, that 
Z l = U^-S- + 0(h s ). The modified output functions are thus y l = Z l (H ), that i 


is, 


v l = u i l *k( H + °^)) + °( ftS ) = y l + D 

Let us now turn to the problem of energy balance under sample and hold. The 
power balance (12.161) implies that 


D a -\-\ H a — 


ah-\-h 


ui(t)y l (t)dt , 


(4.18) 


where we have defined the sampled Hamiltonian H a := H(x a ). A usual way to im¬ 
prove the transient behavior of the system is to add damping by means of a continuous 
control law [23] 


ui(t) = - K ijy 3 (t ) , 


(4.19) 


with {A" ;? } a symmetric and positive semi-definite matrix. With the control law 
(14.191) . the power balance (14.181) results in the dissipation inequality H a+ 1 — H a < 0, 
which guarantees that H a decreases monotonically and, if the right conditions are 
met, the system converges to a state of minimal energy. 
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Suppose that the output is being sampled and that the input is being held at 
intervals of length h. The control sequence is then given by 

u la = ~K ljy i a ( 4 . 20 ) 

and the power balance (14.181) takes the form 

nOth-\-h nh 

H a + i ~H a = u la / y l (t)dt = u lol / y l (ah + r)dr . 

J *h J 0 

Applying Taylor’s theorem to the integral term gives 

H a+1 -H a =J2 «ia (»U + 0(h 2 )) = -hK ljV l a y\ + £ u la O{h 2 ) , 

i i 

so H a decreases when h is small enough and the norm of y a is large enough. 

Since the approximate sampled-data model (14.51) is also PH (cf. item (iii) of 
Theorem 14.2D . it satisfies (again, after applying Taylor’s theorem) 

H a+1 -H a = J2 (y l « h + °( ft2 )) ( 4 - 21 ) 

i 


for some H and y. According to (14.6)) , the energy balance (|4.21|) takes the form 
H a+1 -H a + 0(h s ) = Y, n la (( y l a + 0(h s )) h + 0(h 2 )) . 


The same control sequence (14.201) produces 

H a+1 ^H a = - hK tj y l a y\ + ^ u la O{h 2 ) + 0{h s ) . 

i 

Thus, for s = 2, the qualitative behavior of the approximated sampled data model is 
the same as the exact one: H a decreases when h is small enough and the norm of y a 
is large enough. 

Example: A double planar pendulum (continued). Let us compute a 
sampled-data model for the double pendulum described in the previous examples. 
The first step is to compute a sample-data model for the simple unconstrained PH 
system X'h, u = Dh + uiU i 1 jp-, where H is given by (12.201) and by (12.221) . The 
unconstrained and unactuatecf Hamiltonian vector field Dh describes a pair of masses 
with initial positions r a 0 and r b 0 and initial momenta p a0 p b0 , simply falling under 
the influence of gravity. The exact flow generated by Dh, the drift, denoted by 
(r a +i,p a+ i) = 4>H,h( r <x,Pa), is then given by 


Ck! —|— 1 


a+1 


m a 

h 


Pa x 


„ a h _h 2 

a+1 = r \ + —Pa v a. - 9 — 


-Pb x 


Ck+l 


m a 

h _h 2 

-Pbyoc-9 y 


and 


Pa a;Q!+l Pa x a. 5 


Pa y oc+l = Pay* - ™agh 
Pby *-\-1 Pby* 'W'bgh • 


Pb x *-\-1 Pb x * i 
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Value 


Description 


l a = 0.6[m] 

lb = 0.3[m] 
m a = 0.2 [kg] 
m b = 0.6 [kg] 
g = 9.81[m/s] 


9.81 [m/s] 2 Acceleration due to gravity 


Length of the first link 
Length of the second link 
Value of the first mass 
Value of the second mass 


Table 4.1 

Parameters for the double pendulum. 


The exact flow generated by the control vector field without drift, is denoted 

by (rQ+i,Pa+i) = 4>u,h{r a ,p a )- It is given by 


r Q + l r a 1 * G { a X > 0,y, b x , by } 


and 



a 


a 



b 


b 


From we know that a simple symmetric method of order two for Xjj,u is 


(4.22) 


1pH,u,h 


4*H , if ° fiujh ° 4 > _ 


Notice that <j> u ,h = Id when u = 0, so ipH,u,h = 4 >h ^ ° 4>h ^ = 4>H,h, which is a 
symplectic map because it is the exact solution of a Hamiltonian system. Therefore, 
Xh,u is PH, and, from Theorem 14.21 it follows that the implicit method (14.51) . with 
^H,u,h as in (14.221) is the exact solution of a PH system AH,u,g with Hamiltonian 
function H = H + 0(h 2 ) and output function y = y + 0(h 2 ) (i.e., s = s = 2). 

The sampled-data model was tested using the parameters shown in Table 14.11 
For illustration purposes, we chose a damping control u a = —0.3 -y a and simulated 
the closed-loop system using the sampled-data model (14.51) . Figure lT2l shows the 
discrete-time series of H for different values of h. It can be seen that the time series 
converge and, as expected, the value of H a decreases monotonically when h is small 
enough (in this case, less or equal to 30 [ms]). For comparison purposes, we have 
included the evolution of H that is obtained by simulating (with Matlab’s Simulink) 
the explicit model developed in [4] in series with a sampler and a zero-order hold. 

5. Conclusions. We have extended the second-order integration method pre¬ 
sented in m ■ The original method applies to autonomous Hamiltonian systems and, 
being symplectic, preserves the Hamiltonian structure of the continuous-time system. 
The extended method can be applied to port-Hamiltonian systems, which are Hamil¬ 
tonian systems equipped with input-output pairs. The extended method preserves 
the port-Hamiltonian structure. Affinity in the controls is lost by the method but, 
fortunately, the passivity properties of the continuous-time system can be recovered 
by a suitable redefinition of the output. Interestingly, the relation between the original 
and the new output is also of second order. 

The integration method can be used with the purposes of numerical simulation 
or with the purpose of deriving discrete-time models to be used in the design of 
discrete-time control laws. 
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Fig. 4.2. Results of the numerical experiment. The Hamiltonian function is plotted against 
time. The explicit model (simulated using Matlab’s module Simulink) is compared with the implicit 
model (simulated using a Matlab script). As expected, H is monotonically decreasing when h is 
small enough and the time series converges as h goes to zero. 
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